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Associative learning of temporally disparate events is of fundamental importance for 
perceptual and cognitive functions. Previous studies of the neural mechanisms of such 
association have been mainly focused on individual neurons or synapses, often with an 
assumption that there is persistent neural firing activity that decays slowly. However, 
experimental evidence supporting such firing activity for associative learning is still 
inconclusive. Here we present a novel, alternative account of associative learning in the 
context of classical conditioning, demonstrating that it is an emergent property of a 
spatially extended, spiking neural circuit with spike-timing dependent plasticity and short 
term synaptic depression. We show that both the conditioned and unconditioned stimuli 
can be represented by spike sequences which are produced by wave patterns propagating 
through the network, and that the interactions of these sequences are timing-dependent. 
After training, the occurrence of the sequence encoding the conditioned stimulus (CS) 
naturally regenerates that encoding the unconditioned stimulus (US), therefore resulting 
in association between them. Such associative learning based on interactions of spike 
sequences can happen even when the timescale of their separation is significantly larger 
than that of individual neurons. In particular, our network model is able to account for the 
temporal contiguity property of classical conditioning, as observed in behavioral studies. 
We further show that this emergent associative learning in our network model is quite 
robust to noise perturbations. Our results therefore demonstrate that associative learning 
of temporally disparate events can happen in a distributed way at the level of neural 
circuits. 
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INTRODUCTION 

Associating sequential events happening at different time 
moments is of fundamental importance for a host of perceptual 
and cognitive functions (Wallenstein et al, 1998; Fuster et al, 
2000). One of the important paradigms of such association is clas- 
sical conditioning, in which the pairing of two subsequent stimuli 
is learned such that the presentation of the CS is taken as a predic- 
tor of the subsequent US (Rescorla, 1988). Despite the fact that all 
classical conditioning experiments have an important temporal 
component (i.e., the US occurs after the CS), many theories such 
as the Rescorla- Wagner theory of conditioning do not include 
this timing relationship in their models (Rescorla and Wagner, 
1972; Markram et al, 2011). Furthermore, most of the theo- 
retical models proposed for classical conditioning have focused 
on individual neurons or synapses by assuming the presence of 
slowly decaying firing activity of neurons. This, when combined 
with some synaptic learning rules such as spike timing dependent 
plasticity (STDP), enables associative learning between tempo- 
rally separated events to happen (Gluck and Thompson, 1987; 
Worgotter and Porr, 2005; Drew and Abbott, 2006). However, 



experimental evidence for such slowly decaying firing activity is 
still inconclusive (Ito et al., 2008). 

Here, we present a novel, alternative account of associative 
learning of spatially and temporally separated events, which 
utilizes the emergent dynamics of neural circuits rather than 
individual neurons. Our proposal relies on two basic, known 
neurophysiological features: ( 1 ) spike-timing dependent plastic- 
ity (Abbott and Nelson, 2000; Dan and Poo, 2006), and (2) 
short term depression (STD) (Zucker and Regehr, 2002). In our 
proposal, both the CS and US are encoded as spike sequences 
formed by propagating spiking waves, which emerge from spa- 
tially extended, spiking neural circuits, in which each neuron is 
coupled to nearby neurons in two dimensions, with coupling 
strengths being a function of distance. Associating the CS and 
the US with a time separation greater than 100 ms naturally 
occurs through the timing-dependent interacting dynamics of 
these spiking waves. We further show that our model of asso- 
ciative learning is able to account for the temporal contiguity 
property of classical conditioning, i.e., that the success of asso- 
ciative learning of conditioning is a non-monotonic function 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



July 2014 I Volume 8 | Article 79 | 1 



Palmer and Gong 



Wave interaction based associative learning 



of the amount of time by which the CS precedes the US; in 
particular when the precedence is either too short or too long, 
conditioning is relatively poor, as found in behavioral studies 
(Rescorla, 1988). We also show that our model is robust against 
perturbations. 

The idea of using interacting waves to account for the asso- 
ciative learning of different stimuli can be traced back to Beurle's 
landmark study (Beurle, 1956). Since this work, however, studies 
of association have focused on associatively recalling a stimulus 
based on partial clues, mainly dominated by the view of fixed 
point attractors (Hopfield, 1982). Recently, evidence of propa- 
gating waves in neural systems has been rapidly accumulating 
(Rubino et al, 2006; Benucci et al., 2007; Ferezou et al, 2007; 
Han et al, 2008; Wu et al., 2008; Lubenov and Siapas, 2009; Sato 
et al., 2012), strongly suggesting that it is time to consider the 
fundamental question regarding the potential functional roles of 
propagating waves in the brain (Gong and Van Leeuwen, 2009). 
Our new model of associative learning is a step forward along this 
direction. 

MATERIALS AND METHODS 
NEURAL CIRCUIT MODEL 

We consider a n-by-n two dimensional lattice of integrate-and- 
fire (IF) neurons. We denote the membrane potential of a neuron 
at integer coordinates (i, j) at time t by Vy{t), which has dynamics 
governed by the following equation. 



V,j (t+At) 



e-^'/'V,j it) +y''^+V™(f) ifV(t) < Vth, 
V.j(t)-Vth iiV(t)>Vrh, 



where r = 20 ms is the neural time constant, Af = 1 ms is the 
time step, V™ is a constant external input, V'^"(t) is the input 
received from other neurons at time t. Here both the constant 
external input and the input due to the spiking neurons both have 
an implicit dependence on Af. However, as Af is constant, this 
dependence is absorbed into the appropriate term. A spike is gen- 
erated whenever the neuron reaches a threshold voltage, Vf/, = 1. 
The coupling strength between any two neurons located at (i, j) 
and (/,/) respectively is denoted as Wy./y constructed from a 
"Mexican-hat" function: 



w: 



0 



iidij^^p > Di 



(2) 



where dij, i'f is the Euclidean distance between neurons on 
a lattice with periodic boundary conditions, Ce = 0.4, C/ = 
0.1, dE = VTi, di = ^/42 are constants that determine the shape 
of the coupling function. As interactions between neurons in 
real neural systems have finite ranges, connections between neu- 
rons in the model are constrained to dy,;'/ < -Di = 15. For 
these parameter choices, coupling is excitatory if dy,,'/ < Do = 
y/21\og (4) and inhibitory if Di > dij ff > Do, therefore resulting 
in a lateral inhibitory coupling structure; such lateral inhibi- 
tion could be achieved through disynaptic pathways (Melchitzky 
et al, 2001) or large basket cells (Markram et al., 2004). 
Such spiking neural circuits with lateral inhibition have been 



commonly used to model neural systems including the hip- 
pocampus (Samsonovich and Mcnaughton, 1997). Note that we 
have found that the collective dynamics such as the formation of 
spiking waves in the neural circuit are not sensitive to the values 
of any of the above parameters, as shown in Palmer and Gong 
(2013). 

To enable the total excitatory and inhibitory coupling 
strengths to be changed without varying their spatial extent, 
they are normalized separately using the following normalization 
constants: 

^l = llK''j' if^y,.7^0o, 

^1 = 12 ^'j.>'r if-^o < ^y-'7 < -^1- (3) 

,'/ 

The coupling strengths are then given by 



'i.'f 
wf. 

'J.'/ 



iidij^i'jt < Do, 

if Do < dy_,y < Di, 
if%,7 > Di, 



(4) 



where We =1.6 and IV; = 2.1 are the total excitatory and 
inhibitory coupling strengths, respectively. 

All synaptic connections between neurons are subject to both 
long-term and short-term plasticity; for long-term plasticity, we 
incorporate STDP for all synapses as follows: 



AWy,,y (f) 



tl<t t'.,.,<t 



(5) 



ii djj fji > Di , 



where f^ is the kth firing time of neuron (;,;) and t'-,j, is the zth 
firing time of neuron (i' , j'), the sum is over all spike pairs before 
time f, and H ^f^, f|y^ is the STDP window function: 



H(fj 



■post ) '■pre j 



A+ exp (-Af/T+) if Af > 0, 
-A- exp (Af/T_) if Af < 0, 
0 ifAf = 0, 



(6) 



where fpre and tpost are the firing times of the pre- and post- 
synaptic neurons, respectively, and Af = tpost — tpre- The STDP 
timescales, r_|- = t_ = 20 ms, are similar to values found in 
experimental studies (Bi and Poo, 1998; Zhang et al., 1998), 
and A+ = A_ = 0.00025, making H an odd function in Af. We 
incorporate the STDP into the coupling strength to obtain the 
time-dependent coupling strength (Palmer and Gong, 2013): 



(f) = Wy,,y + AWy,,y(f). 



(7) 



STDP may cause some coupling strengths to reach values that are 
not within biological limits; to prevent this, STDP requires a hard 
limit of synaptic modification (Song et al, 2000; Izhikevich et al.. 
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2004). No further constraints are required. We limit changes to 
coupling strengths as follows: 

|(1 - AW„^)Wy,,y (0)| < I (f)| 

< |(l+AW„,„^)Wy,,y(0)|, (8) 

where Wy^iy (0) is equal to the initial connection strengths as 
described in Equation (4), and Wmax = 0.18 controls the max- 
imum possible amount of STDP. The results discussed later are 
not sensitive to this value. 

Short term plasticity, also called dynamical synapses, refers to 
a phenomenon in which synaptic efficacy changes over time such 
that it reflects the spiking history of a neuron. In contrast with 
long-term plasticity dynamics such as STDP, short term plastic- 
ity induces temporary modification to synaptic efficacy; that is, 
without continued presynaptic activity, the synaptic efficacy wiU 
return to its baseline level. We incorporate such short term plas- 
ticity by using the model described in Markram et al. (1998); in 
this model the synaptic efficacy A^,j, of the neuron at {i', /) is 

4f = (9) 

4f = ^+«?7 ' ( 1 - C/) exp (- Af,- ) , (10) 

r^, = l + (r^-i-r^i„^-i-l)exp(-A*;-VD), (H) 
u],j, = U,rl.,= l (12) 

where Ajj^ represents the absolute synaptic efficacy of the fcth spike 
of the neuron, u^,-, represents the utilization of synaptic efficiency, 

represents the availability of synaptic efficiency, A^,^ ^ is the 
time between the k — 1th and kth spike of the neuron located at 
(f,/) and -F = 0.005 s, D = 0.11 s are the timescales represent- 
ing the recovery of the u and r variables, respectively (Markram 
et al., 1998). Since D JS> _F, these synapses are depressed, i.e., 
the amplitudes of excitatory postsynaptic potentials caused by 
closely-spaced successive spikes decrease over time; this is referred 
to as short term depression (STD). Furthermore, the closer in 
time successive spikes are, the smaller the efficacy of the spikes, 
with synaptic efficacy recovering over a timescale which is ~ D. 
Combining STDP described in Equation 7 with STD described in 
Equations 9-12, the input the neuron (z,;) receives from other 
neurons to which it is coupled is: 

Vlj'(t) = J2w,jjj,it)A'l,^,Q'){V,,j,-V,h), (13) 
"7 

where 0 {x)is the Heaviside step function: © (x) = 1 if x > 0, and 
© (x) = 0 if X < 0. 

Since real neural systems are full of noise, we also add noise 
to the model. We incorporate noise through adding noisy spon- 
taneous firing; for each neuron, a spike is emitted with small 
probability to generate random sparse spontaneous firing with 
a frequency around 0.1 spikes per neuron per second, which 
is within the range of sparse spontaneous firing rate as found 



in cortical neurons (Griffith and Horn, 1966; Koch and Fuster, 
1989). 

RESULTS 

ASSOCIATIVE LEARNING BASED ON THE INTERACTIONS OF SPIKE 
WAVES 

To model the association of the conditioned and unconditioned 
stimuli, we first demonstrate how each stimulus can be repre- 
sented in the network. To this end, we add a brief localized 
stimulus representing either the CS or the US to neurons at 
different locations in the network. These stimuli are able to 
evoke spiking wave patterns as shown in Figure lA. As the wave 
pattern propagates across the network, it naturally gives rise to 
a spike timing sequence; each stimulus, therefore, is encoded by 
a spike sequence. In the present study, we use the terms spike 
sequence and propagating spiking wave interchangeably. This 
method of generating spike sequences utilizing propagating pat- 
terns in recurrent neural circuits is similar to that previously used 
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FIGURE 1 I Propagating spiking waves in tlie neural circuit model can 
be used to implement associative learning. (A) A wave propagating 
through the neural circuit with V" = 0.0429. Color shows the value of 
Vjjit). The wavefront and its tail with reduced potential are clearly visible. 
(B) Schematic of wave interactions for associative learning. Black boxes 
show the location of the stimuli where the CS is initialized at the top of the 
network and the US on the left. Gray regions indicate the paths of the 
propagating spiking waves, the red region shows the interaction region, 
and the blue region shows the response area. (C) Schematic of wave 
regeneration. The wave encoding the CS (red wavefronts, black path) is 
initialized at the top of the network and begins propagating toward the 
interaction region. When it enters the interaction region, it regenerates the 
US (green wavefront, blue path), which was initialized at the left of the 
network, thereby achieving associative learning. (D) Number of successful 
associations in 60 trials. There is a sharp increase in the number of 
successful associations with a peak near 100 ms, followed by a gradual 
decline. To ensure that the predictive relationship is maintained, we test to 
see if the US can regenerate the propagating wave encoding the CS. 
Bolded points show when the number of responses in this case is 
significantly lower (p < 0.01 using a f-test) than the number of successful 
associations; this indicates that the predictive relationship between the CS 
and the US has been maintained. 
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in Itskov et al. (201 1), in which spike sequences formed by local- 
ized moving patterns were used to encode elapsed time in behav- 
ioraUy relevant contexts. Representing external stimuli utilizing 
spike sequences in our model is consistent with a growing body 
of experimental results showing that such sequences pay an 
important role in a host of perceptual and cognitive processes 
(Pastalkova et al, 2008; Bathellier et al., 2012; Harvey et al, 2012; 
XuetaL, 2012). 

As shown in Figure IB, both the CS and the US can evoke 
propagating waves which naturally produce spike sequences; in 
this example, the CS is represented by the spiking wave traveling 
from top to bottom, and the US is represented by the spiking wave 
traveling from left to right. Here we assume that when the wave 
evoked by the US reaches the corresponding blue box on the right 
of Figure IB, a response corresponding to the US is generated; 
we do not explicitly model how the direction of wave propaga- 
tion can be read out to generate a behavioral response. However, 
such read out could be achieved through spatial arrangement of 
dendritic receptor fields. Indeed, it has been found that dendrites 
are sensitive to spike sequences arriving from different directions 
(Branco et al, 2010); in a recent modeling study, this property 
has been used to read out information contained in propagating 
waves (Heitmann et al., 2013). 

We now demonstrate that associative learning of the CS and 
the US naturally arises in the network with STDP and STD. 
Although there has been a significant number of experimental 
studies measuring how such learning can occur in a variety of 
situations, the general protocol that has been used is: the two 
paired stimuli are presented repeatedly to allow learning to hap- 
pen; the CS is presented first. After a delay time, Af, which is 
constant for all learning trials, the US is then presented. After 
training, when the CS is presented, the response of the US is gen- 
erated, indicating that the CS can be taken as a predictor for the 
US (Rescorla, 1988). Using the same protocol as used in behav- 
ioral studies, we apply the paired stimuli to the network model 
for 60 trials. After some learning trials (typically < 20), when the 
CS is presented alone, it is able to generate the unconditioned 
response. As illustrated in Figure IC, this occurs because when 
the spike sequence (wave) encoding the CS propagates along 
its path, it is able to regenerate the spike sequence (wave) that 
encodes the US. 

THE TEMPORAL CONTIGUITY PROPERTY OF CONDITIONED 
ASSOCIATIVE LEARNING 

As shown in behavioral studies, the success of associative learning 
of conditioning is a non-monotonic function of the amount of 
time. At, by which the CS precedes the US (Rescorla, 1988). When 
Af is either too short or too long, conditioning is relatively poor. 
However, there is some intermediate set of time intervals when 
conditioning is most successful. This temporal contiguity prop- 
erty of classical conditioning has been found to be invariant across 
many different experiments covering a wide variety of stimuli, 
although there is no absolute temporal interval that produces 
the best conditioning. For instance, a similar non-monotonic 
functional form has been found when testing both rabbit eyelid 
response and pigeon keypeck response, however, the time inter- 
vals for optimal conditioning were on the order of 100 ms and 5 s. 



respectively (Smith et al., 1969; Barker and Smith, 1974; Rescorla, 
1988). 

We now demonstrate that this temporal contiguity property 
can be reproduced by our propagating wave-based model of asso- 
ciative learning. When both STDP and STD are included, we run 
60 trials of the paired US and CS. As shown in Figure IC, for 
some values of Af, after training it is possible for the CS to gen- 
erate the response of the US. In our testing regime, we apply the 
CS alone after each learning trial. If the US response is generated, 
it is counted as a successful association. As in behavioral stud- 
ies (Rescorla, 1988), we count the number of successful trials, and 
repeat this procedure for different A f values to study how success- 
ful conditioning varies with Af. Figure ID shows that our model 
is able to reproduce the temporal contiguity property of condi- 
tioning, as observed in experimental studies (Smith et al., 1969; 
Barker and Smith, 1974; Figure 1 of Rescorla, 1988). In particu- 
lar, the model shows a very low number of responses for small Af 
followed by a sharp increase with a peak at Af ^ 100 ms, which 
is then followed by a decay in number of responses for larger Af. 
The associative learning therefore is relatively poor if A f is either 
too short or too long, and it appears to be a non-monotonic 
function of Af. Note that despite the general trend of the decay 
of successful learning for larger Af (the mechanism underlying 
such decay is illustrated below), a significant number of success- 
ful learning can stiU happen even when Af ~ 1 s, a time scale that 
is much larger than that of individual neurons or synapses. 

MECHANISMS OF ASSOCIATIVE LEARNING BASED ON SPIKING WAVE 
INTERACTIONS 

As demonstrated above, associative learning of the CS and the 
US can naturally occur in our model, utilizing spike sequences 
generated in the network. The question that arises naturally is: 
what are the mechanisms underlying such associative learning? 
It has been proposed that spike sequences are important neural 
substrates for perceptual and cognitive functions (Abeles, 1991; 
Kumar et al., 2010). However, the timing relationship between 
such sequences has only been explored for their zero-lagged 
synchrony and its potential computational roles for feature bind- 
ing (Abeles et al, 2004). As proposed recently, cell assemblies 
organized as spike sequences could be the basic tokens of the 
"neural syntax" (Buzsaki, 2010), just as words are the basic tokens 
of the syntax of language; clearly, such complex neural syntax 
requires that different spike sequences exhibit timing relation- 
ships more complicated than those of zero-lagged synchrony. We 
now demonstrate that the spike sequences formed in our neu- 
ral circuits have complex interactions that are timing dependent 
and that such timing-dependent interactions of spike sequences 
are essential for implementing associative learning at the level of 
neural circuits. 

To study how the interactions between the propagating spik- 
ing waves evoked by the CS and the US depend on their time 
separation, Af, we systematically vary Af from 0 to 1000 ms. 
Since the CS is added to the network prior to the US, for clar- 
ity we label the spiking waves evoked by the CS and US as the 
first and second waves, respectively. We first investigate the behav- 
ior of the network in the first learning trial as AW^ ,-y (f) in 
Equation 5 is small, and the additional dynamics added by STDP 
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FIGURE 2 I Interacting dynamics of spiking waves. The propagating path 
of the wave encoding the CS is denoted by the red line; that of the wave 
encoding the US is denoted by the black line. (A) A repulsive interaction 
(Af = 0 ms). (B) A suppressive interaction (At = 20 ms). (C) The second 
wave passes through the tail of the first (At = 360 ms). 



are minimal, which simplifies analysis. As shown in Figure 2, the 
interactions can be largely grouped into three main categories, 
namely, repulsion, suppression and pass-through. The repulsive 
interaction occurs when the second wave interacts with either the 
front or the tail of the first wave, causing the paths of the both 
spike waves to be altered; this typically occurs for small values of 
Af (Af < 15 ms). For larger values of Af (15 < Af < 100 ms), 
suppression can occur when the second wave interacts with the 
tail of the first one, as a result the second wave is annihilated. 
For both the repulsive and suppressive interactions, it is appar- 
ent that the wave evoked by the US is unable to pass through the 
path of the one evoked by the CS; i.e., the propagation path of 
the first wave is blocked by the tail of the second wave. However, 
when Af > 100 ms, the second wave is able to pass through the 
path of the first one. These wave interactions therefore depend 
on their temporal separation; in other words, their interaction 
dynamics are timing dependent. Similar timing-dependent inter- 
action dynamics of evoked propagating waves has been found 
in the rat visual cortex; as reported in Gao et al. (2012), when 
the inter- stimulus interval was between 80 and 300 ms, a second 
propagating wave was significantly suppressed by the first one. 
This is similar to the suppressive interactions shown in our model 
and there is certain overlap between the interval (80-300 ms) 
found in Gao et al. (2012) and that in our model (15-100 ms). 
For inter- stimulus intervals greater than 300 ms, the amplitudes 
of the two waves showed no significant difference (Gao et al., 
2012). This is similar to the pass-through that we observed for 
large values of Af, although for our model this behavior occurs 
when Af > 100 ms. 

THE ROLES OF STOP AND STD IN ASSOCIATIVE LEARNING 

In the previous Section, we have illustrated the dynamics of 
interactions between multiple propagating spiking waves when 



AWy ;y (f) is small (i.e., the beginning of the learning process). 
After allowing learning to proceed, when Af is small, the block- 
ing effect still exists as shown in Figures 2A,B. When Af is large, 
the second wave can freely pass through the tail of the first one 
(Figure 2C); however, after several trials of learning, the presence 
of the first wave along can regenerate the second one, as shown 
in Figure IC. We now show that this occurs because STDP plays 
a crucial role in learning the paths of the propagating waves and 
their interactions. 

We first examine how coupling strengths are changed as a 
result of STDP during the learning process of the paired CS 
and US. As a wave pattern propagates across the neural circuit, 
the coupling strengths aligned with the direction of its propa- 
gation will be increased as STDP increases coupling strengths 
when presynaptic neurons fire before postsynaptic ones. On the 
other hand, as the STDP window function (Equation 6) is odd, 
connections in the direction opposite to wave propagation will 
be decreased, as shown in Figure 3A. Whilst coupling strength 
changes resulting from a wave pattern traveling over a path once 
are relatively small, if it repeatedly travels over the same path (as 
in the training process) significant changes to coupling strengths 
can occur, with the maximum change limited by the bounds 
to STDP (Equation 8). As shown in Figure 3A, it is simple to 
understand coupling strength changes due to a single pattern, but 
analyzing the changes to connection strengths within the inter- 
action region are more complex, as each neuron has around 700 
connections and multiple wave patterns travel over each neuron 
from different directions. To simplify, we define a function that 
can be used to quantify changes to the strength of synaptic con- 
nections from a single neuron as a function of a single variable, 
namely the angle of the connection on the two-dimensional, spa- 
tially extended spiking neural circuit. This is achieved by creating 
a new polar coordinate system centered on the neuron of interest 
and averaging connection strengths over certain angular segments 
(Figures 3A,B). The change of coupling strength as a function of 
angle is: 



A,-; (0) = max 



@{d- 



.DlA%>+xy+;(f) 



tan 



-'(1)1) 



V 



e( 



@(d-\e 



— tan" 



'(f)l) 



(14) 



where d is a small number that controls the angular resolution 
of Ay (9), Di is the maximum coupling distance (Equation 2), 
and X and y index over aU neurons that the neuron at location 
is connected to. Ay (6*), defined in Equation 14, measures 
the amount of synaptic changes caused by STDP as a function of 
angle, 6, from the presynaptic neuron (Figure 3A). This is cal- 
culated by summing over all connected neurons and then using 
a Heaviside step function to ignore those neurons which are con- 
nected to the main neuron at angles different from 0, i.e., neurons 
outside the red sector in Figure 3A. As neurons are aligned along 
a regular grid, the spatial density of connections is uniform but 
the angular density of connections is non-uniform; to correct for 
this, we normalize Ay (6) by dividing it by the number of neurons 
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FIGURE 3 I Changes to synaptic coupling strengths due to STDR (A) 

Changes in coupling strength (AlA/j- /'/(f)) after a wave pattern has traveled 
over a neuron. Coupling strengths in the direction of propagation are 
increased, and strengths in the opposite direction are decreased. To 
calculate A,y(9), the coupling strength changes within the circular sector 
(red lines) centered at an angle of 0 (black lines) are averaged to give the 
value of Ajj{9). (B) Aj- (9) for the neuron in (A). There is a clear peak pointing 
in the direction of wave propagation across the neuron. (C) A,; (6) for the 
20 X 20 sub-grid of neurons near the interaction region. The original paths 
of the two wave patterns corresponding to the CS and the US are marked 
with black and red lines, respectively. The gray-filled circles mark neurons in 
the interaction region. The polar plots are laid out so that each plot 
corresponds directly to a single neuron in the grid of neurons. For scale, the 
black circles show a mean change in connection strength of 10"". In the left 
and top of the figure, only a single wave passes over each neuron, leading 
to connection increases in only one direction, similar to (B). The areas in 
the corners with no connection strength changes because in the regions 
the wave does not pass over In the interaction region, neurons typically 
show an increase in coupling strength in two directions, corresponding to 
the top-to-bottom CS and the left-to-right US. The resulting bimodal 
distribution of coupling strengths enables the spontaneous wave 
regeneration to happen. The red-filled circle shows a neuron on the edge of 
the path. As A,y(e) is slightly biased toward the center of the wave pattern, 
waves initialized near this neuron will quickly converge to the learned path. 



that we are averaging over. The negative component of A Wj, ,y is 
removed for simplicity, as it is almost always in the exact reverse 
direction when compared to the positive component. 

As discussed above, we expect that when a single wave passes 
over a neuron, A,j(d) will have a single peak that is approxi- 
mately aligned with the direction of propagation of the wave. This 
is shown in Figure 3B where Ay(S) shows a single peak, which 
is aligned with the propagation direction of the wave. We now 
calculate A/,- (6*) for all neurons near the interaction region after 
training by repeatedly applying the paired US and CS stimuli, as 



shown in Figure 3C. We find that in this case, the behavior of 
Aij{0) near the interaction region is more complex than the case 
when only a single wave is present (Figure 3B). In particular, as 
waves are interacting within this region, both wave patterns cause 
changes to coupling strengths in their directions of propagation. 
As a result, neurons within the interaction region show two peaks 
in Ay (6), corresponding to the directions of the two propagating 
wave patterns (Figure 3C), namely, top-to-bottom and left-to- 
right. This bimodal distribution of coupling strengths developed 
during the learning process endows the network with the prop- 
erty of spiking wave pattern regeneration as shown in Figure IC. 
Specifically, when the wave encoding the CS propagates alone and 
enters the interaction region (Figure 3C), neurons to the right of 
the wave pattern wiU receive significant excitatory input due to the 
increased coupling strengths in that direction (Figure 3C). This 
will cause certain neurons in that direction to fire; if the num- 
ber of these firing neurons reaches a critical value, it is possible 
for a wave pattern propagating along the right-to-left direction to 
be regenerated and this regenerated wave will then travel along 
the path of the wave evoked by the US. This wave regeneration 
allows for the wave evoked by the CS to generate the uncondi- 
tioned response, a key feature in our associative memory model 
(Figure IC). Here it is interesting to note that in our theoretical 
account of associative learning, the neurons in the intersection 
region play an essential role in associating the US and the CS 
together; this notion is similar to the concept of "nodal" neu- 
rons discussed in Eichenbaum et al. (1999), in which it was 
hypothesized that some nodal neurons could encode the com- 
mon features of different temporal events, so that these temporal 
events could be associated to form episodic memory (see Figure 3 
in Eichenbaum et al., 1999). 

Aside from changing coupling strengths during the training 
process (Figure 3), STDP also has a notable secondary effect, 
namely, that the speed of spiking wave propagation undergoes 
an increase of 25% after only 5 training trials. This property is 
qualitatively similar to the experimental data reported in Xu et al. 
(2012), in which visual stimulus (a moving dot) could induce 
spike sequences and the propagation speed of these sequences 
increased after a training process, i.e., many repeated presen- 
tations of the stimulus. Such increased speed of training or 
experience-induced spike sequences has also been observed in the 
hippocampus and in the neocortex (Ji and Wilson, 2006; Euston 
et al., 2007). 

In associative learning of classical conditioning, there is a pre- 
dictive relationship between the CS and the US. This allows the CS 
to predict the US, whilst the US does not predict the CS (Rescorla, 
1988). In our model, we find that with STDP but without STD, 
this predictive relationship is violated when At ^ 200 ms, as the 
second wave encoding the US is capable of generating the CS 
response. This happens because without STD, when the second 
wave approaches the interaction region, the membrane potential 
of neurons along the paths of the first and the second wave has 
been restored to their initial values, thus the membrane potential 
would be near-identical for both the first and the second waves. 
Incorporating STD, however, breaks this symmetry by adding 
an extra slowly-recovering effect, namely synaptic efficacy. As 
shown in Figure 4, the reduction in synaptic efficacy caused by 
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the propagation of the first wave persists long after the neuron 
membrane potentials have recovered to their equilibrium values. 
As described in Equation 13, synaptic efficacy can affect the total 
inputs to neurons; namely, the smaller synaptic efficacy is, the 
smaller the inputs to other neurons are. This means that when 
the second wave reaches the interaction region (Figure 3C), the 
inputs caused by its spiking wave front to the neurons along the 
direction of the first wave are reduced due to the small synaptic 
efficacy caused by the first spiking wave. In this way the sec- 
ond wave is unable to regenerate the first wave; therefore, the 
predictive relationship between the CS and the US is maintained. 

DECAY OF SUCCESSFUL ASSOCIATION DUE TO NOISY SPONTANEOUS 
FIRING ACTIVITY 

The general function form of the temporal contiguity prop- 
erty of classical conditioning is that the success of association 
between the CS and the US is a non-monotonic function of their 
time separation. At, as found in experiments (Rescorla, 1988) 
and reproduced in our modeling study (Figure ID). In particu- 
lar, for large values of Af, the amount of successful association 
decreases as Af increases; this decay of association over time 
is ubiquitous for different experimental protocols (see Figure 1 
of Rescorla, 1988), but the underlying neural mechanism stiU 
remains unclear. We now show that in our model, this decay of 
association occurs due to the random modifications of synap- 
tic coupling strengths caused by spontaneous firing activity (see 
Materials and Methods). 

As demonstrated above, during the training process, the 
synaptic coupling strengths along the propagation paths of the 
first and second waves are modified by STDP to store these spik- 
ing waves. However, due to the random nature of the spontaneous 
firing activity, these coupling strengths wiU be changed in a ran- 
dom way; i.e., the coupling strengths will be randomly potentiated 
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FIGURE 4 I Comparison of the recovery process of the membrane 
potential (V) and synaptic efficacy {A). A propagating spiking wave is 
initialized and the value of l/and A is measured along its path. Both are 
normalized with their resting values equal to 1 enabling them to be plotted 
on the same scale for comparison. Note that the x-axis is a log scale, whilst 
the y-axis is linear. The front of the wave is near x = 10, at the peak of the 
green line which indicates a neuron emitting a spike. It is clear that the 
membrane potential has recovered to close to the initial value within ajlOO 
grid units of the wavefront, whilst the synaptic efficacy takes «;500 grid 
units to recover. 



or depressed through STDP, depending on the temporal order 
of the random spikes happening to two coupled neurons. One 
would then expect that such random changes can make the cou- 
pling strengths to behave like a random walk. To test this, we 
calculate the variance of the random change of coupling strengths 
AWij i'ji(t) as a function of time. Figures shows that that the 
variance increases linearly over time, i.e., (t) oc f, which is a 
characteristic feature of a random walk. As in Fusi and Abbott 
(2007), we then define the memory traces of the propagating 
spiking waves stored in the modified synapses as "signal," S. The 
random walk of coupling strengths due to random spontaneous 
firing activity would cause the memory traces to fluctuate and 
degrade; such fluctuation introduces "noise," N to the memory 
traces, defined as the standard deviation of the random changes to 
the signal. After the training process, the long-term modifications 
of synaptic strengths due to the propagating waves are AWy 
so the signal S oc AWy j-y. Due to its random walk nature, the 
random modification of synaptic strengths N oc t^^^. The quan- 
tity that can be used to characterize the quality of the stored 
memory of propagating waves is the signal-to-noise ratio, S/N, 
and S/N oc AWy ,-y/f^/^. The signal-to-noise ratio thus decays 
monotonically as Af increases, meaning the memory traces of 
propagating waves decay monotonically; this therefore results in 
a reduction in the number of successful associations between the 
two stimuli. 

ROBUSTNESS OF ASSOCIATIVE LEARNING 

We now show that our associative learning implementation is 
robust to external perturbations. The propagating waves used to 
encode the CS and the US have wavefronts that are composed 
of many firing neurons, which make them robust to some per- 
turbations. In particular, when a propagating wave front is only 
partially activated by external perturbations, it can recover to its 
original size; this is illustrated in Figure 6A, where an obstruc- 
tion is created by forcing blocks of neurons to zero potential. 
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FIGURE 5 I Variance of random synaptic modifications over time. This 
variance increases linearly over time, i.e., (t) oc t, indicating AW/,- /y (f) 
undergoes a random walk. For the linear fit, the coefficient of 
determination, = 0.99997. 
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FIGURE 6 I Propagating spiking waves are stable with respect to 
perturbations. (A) The size of traveling wave patterns is stable. We 
initialize a wave pattern at the top of the image. It then propagates toward 
two regions that are fixed with Vij(t) = 0 (the two blue rectangles). After 
traveling past these obstacles, the wave pattern quickly returns to its 
original size. This demonstrates that the size of the wave pattern is stable. 
(B-E) STDP can enable a path to be learned (schematic). (B,C) show initial 
training of a path without and with STDR respectively. (D,E) Show the 
behavior when a test pattern is initialized to the side of the original pattern. 
With STDP (E), the pattern converges to the original path. Without STDP 
(D), the pattern continues in a straight line. 



i.e., Vy (f) = 0. After propagating through the obstacle, the wave 
pattern quickly returns to its original size, before continuing to 
propagate along its original path. The robustness of the size of 
the wave pattern can be understood through the shape of the 
coupling function (Equation 2). We note that the wave pattern 
size is close to the minimum distance for inhibitory connections 
Dq. As a result, if the wave pattern expands, the amount of inhi- 
bition increases, allowing the pattern to return to its original 
size. Similarly, if the pattern shrinks, the amount of inhibition 
decreases, allowing the pattern to return to its original size. 

Furthermore, the model is robust with respect to perturba- 
tions in the CS and US locations. We consider two cases; one with 
STDP and another without STDP to illustrate the importance 
of STDP in correcting such perturbations. We first run training, 
as shown in Figures 6B,C, by presenting a stimulus at the same 
location for both cases. After training, the location of the wave ini- 
tialization site is shifted slightly to a different side to test whether 
it can return to the learned path, as shown in Figures 6D,E. 
When STDP is included in the network (Figure 6E), the test pat- 
tern quickly returns to the original learned path, indicating that 
for this case, the learned spiking wave is robust to the changes 
of the initial stimulus location. Such perturbations in the ini- 
tial locations can be regarded as errors in the initial stimulus. 
The dynamical behavior of returning to the learned path there- 
fore indicates that the dynamics of the wave patterns now have 
an error-correcting property. This error-correcting property is 
analogous to that shown in conventional attractor neural net- 
works (Hopfield, 1982), in which perturbations to state away 
from a learned attractor are corrected as the state quickly returns 
to the attractor. To understand how STDP has enabled such 
error correction, we investigate how STDP has affected the cou- 
pling strengths of neurons as waves propagate across them, as 



shown in Figure 3C; in particular, the red- filled circle in the figure 
represents a single neuron located near the edge of a wave pattern. 
We observe that Ay(f9) is angled toward the center of the wave 
pattern. As a result of these modified coupling strengths, neurons 
on the original path receive extra excitation. This causes the wave 
to return to the learned path, as shown in Figure 6E. However, 
without STDP, the test pattern does not return to the previous 
path; instead, it propagates along a different path (Figure 6D). 

DISCUSSION 

In this study, we have proposed a novel theoretical account of 
associative learning of temporally disparate events in the con- 
text of classical conditioning. In our model, temporally separated 
events such as the CS and the US are encoded by different spike 
sequences formed from propagating spike waves. After training, 
the sequence encoding the CS can regenerate the sequence encod- 
ing the US, therefore enabling successful association between 
them. Associative learning therefore happens in a distributed way 
at the level of neural circuits, without slowly decaying firing activ- 
ity of neurons. Furthermore, as we have demonstrated, our model 
is able to account for the temporal contiguity property of classical 
conditioning and is robust to certain types of noise perturbations. 

In our model of associative learning, the timing relationships 
of firing activity of groups of neurons, including spike sequences 
formed from these neurons and their timing-dependent interac- 
tions are essential for associative learning. Spike sequences have 
indeed been found to be important for a range of perceptual 
and cognitive processes (Pastalkova et al., 2008; Bathellier et al, 
2012; Harvey et al, 2012; Xu et al, 2012). However, most of 
these studies have focused on just one such sequence and higher- 
level timing relationships between different sequences have 
not been explored. Similarly, previous modeling studied have 
mainly focused on modeling the formation of spike sequences 
(Abeles, 1991; Kumar et al, 2010); interactions between different 
sequences have been rarely discussed and such discussion has just 
focused on their zero-lag synchrony and its role in feature binding 
(Abeles et al, 2004). However, rich timing-dependent interac- 
tions of different spike sequences, as illustrated in our model 
in the context of associative learning, could be of fundamental 
importance for brain functions. 

In our model of associative learning of classical conditioning, 
propagation of spiking waves across neural circuits is the mech- 
anism underlying the formation of spike sequences. As we have 
illustrated, such propagation of neural activity is essential to make 
local information available at larger spatial and temporal scales, 
therefore associating signals that are distributed over space and 
time. Interestingly, the notion that neural activity propagation 
is important for classical conditioning is reminiscent of Pavlov's 
original speculation that "automatic irradiation" of neural firing 
activity evoked by external stimuli could be the neural mecha- 
nism underlying classical conditioning (Pavlov, 1927). Recently, 
evidence of propagating waves in the brain has been accumu- 
lating very rapidly (Rubino et al, 2006; Benucci et al, 2007; 
Ferezou et al., 2007; Wu et al, 2008; Lubenov and Siapas, 2009; 
Sato et al., 2012). In Ferezou et al. (2007), it was speculated that 
multiple waves activated in a distributed manner may be essen- 
tial for "associative plasticity and learning." In Gong and Van 
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Leeuwen (2009), it was proposed that propagating wave patterns 
are basic building blocks for neural circuits to carry out dis- 
tributed dynamic computation; in this paradigm, information 
is encoded in these dynamic patterns and communicated to the 
different parts of the cortex based on their propagation, and infor- 
mation is processed when they collide or interact with each other. 
Here, we have presented such a model for associative learning 
based on interacting spiking wave patterns. 

Interacting waves have been used to account for classical con- 
ditioning in Beurle (1956), in which it was shown that the CS 
and the US can evoke two plane waves. After learning, which was 
implemented through modulating firing thresholds of neurons, 
the activation of one wave is able to regenerate another wave, 
referred to as "pulse regeneration" (Beurle, 1956). In this way, an 
association between two events can be formed. In Beurle's study, 
two waves which occurred at the same time were used to model 
association between the CS and the US; this is in contrast to real 
psychophysical data showing that when the CS and the US occur 
together, no association can be formed (Rescorla, 1988). In our 
study, however, the temporal separation between the US and the 
CS is explicitly modeled and a salient feature of classical condi- 
tioning, namely the temporal contiguity of classical conditioning 
can be reproduced. Furthermore, in our model, the learning of 
wave propagation paths is achieved by using synaptic dynamics 
such as STDP and STD that have been widely observed in many 
different parts of the brain (Abbott and Nelson, 2000; Chung 
et al., 2002); in Beurle's original work, however, there was no 
explicit synaptic plasticity considered and experimental evidence 
of learning achieved by modulating firing threshold is rare. 

Both STDP and STD play important roles in our model of 
associative learning. As we have illustrated, STDP allows spike 
sequences to be learned and importantly such learned sequences 
are very robust against different types of perturbations, including 
changing the location of the initial stimuli. This robust learning of 
spike sequences can occur because asymmetric synaptic weights 
developed in our model enable perturbed waves to converge back 
to a nearby learned path. As we have illustrated, these asymmet- 
ric weights develop due to the temporal asymmetry of STDP and 
because waves move sequentially across the circuit. STDP also 
pays an important role in associating different events together 
by increasing the coupling strengths along the propagating paths 
of different spiking waves to form an intersection point between 
them. The role of neurons around the interaction points estab- 
lished by STDP in our model of associative learning is similar that 
of "nodal cells" introduced in Eichenbaum et al. (1999), in which 
it was proposed that some hippocampal neurons may function 
like nodal cells to link sequential events to form a memory space. 

Short term depression (STD) is another essential element in 
our model of associative learning. Without STD, the US and the 
CS would experience the network in an essentially identical state, 
therefore they are likely to generate similar responses. Previous 
models have broken this symmetry by assuming the existence 
of persistent neural firing activity that decays slowly (Drew and 
Abbott, 2006). However, the existence of such input is still incon- 
clusive (Ito et al, 2008). In contrast, with STD the reduction in 
synaptic efficacy, after the wave corresponding to the CS propa- 
gates through the network, causes the symmetry between the US 



and the CS to be broken. In this model of associative learning, 
therefore, there is no need to assume any persistent neural firing 
activity. Since generating spikes is energetically expensive for neu- 
rons (Niven and Laughlin, 2008), our model potentially provides 
an economical and robust way to implement associative learning. 

Since our model for associative learning is based on several 
general properties of neural circuits, including the existence of 
stimulus evoked propagating waves, STDP and STD, it may be 
applicable to a variety of neural systems. Examples to which 
our model could be applied include the hippocampus, which is 
generally regarded as a place where association of spatially and 
temporally disparate events can form (Rawlins and Olton, 1982; 
Eichenbaum et al., 1999). Notably, both spike sequences and 
propagating waves such as propagating waves of theta oscillations 
(Lubenov and Siapas, 2009; Patel et al., 2012) and propagating 
waves of sharp wave ripples have been observed in the hippocam- 
pus (Patel et al, 2013). Another area in neural systems, to which 
our model can be applied is primary visual cortex, in which it 
has recently been found that moving stimuli can evoke spike 
sequences and these sequences can be recalled once a clue is pre- 
sented (Xu et al, 2012). In these neural systems, it would be 
interesting to investigate timing-dependent interactions between 
different spike sequences and their functional roles. 
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